Differential thalamocortical interactions in slow and fast spindle generation: A computational model

Cortical slow oscillations (SOs) and thalamocortical sleep spindles are two prominent EEG rhythms of slow wave sleep. These EEG rhythms play an essential role in memory consolidation. In humans, sleep spindles are categorized into slow spindles (8–12 Hz) and fast spindles (12–16 Hz), with different properties. Slow spindles that couple with the up-to-down phase of the SO require more experimental and computational investigation to disclose their origin, functional relevance and most importantly their relation with SOs regarding memory consolidation. To examine slow spindles, we propose a biophysical thalamocortical model with two independent thalamic networks (one for slow and the other for fast spindles). Our modeling results show that fast spindles lead to faster cortical cell firing, and subsequently increase the amplitude of the cortical local field potential (LFP) during the SO down-to-up phase. Slow spindles also facilitate cortical cell firing, but the response is slower, thereby increasing the cortical LFP amplitude later, at the SO up-to-down phase of the SO cycle. Neither the SO rhythm nor the duration of the SO down state is affected by slow spindle activity. Furthermore, at a more hyperpolarized membrane potential level of fast thalamic subnetwork cells, the activity of fast spindles decreases, while the slow spindles activity increases. Together, our model results suggest that slow spindles may facilitate the initiation of the following SO cycle, without however affecting expression of the SO Up and Down states.


Introduction
Sleep plays an important role in memory consolidation. Two thalamo-cortical oscillatory rhythms, the sleep slow oscillations (SO) and sleep spindles of non-rapid eye movement (NREM) sleep play a comprehensive role in declarative memory consolidation [1][2][3]. It is believed that cortical SOs provide a temporal frame where recently-acquired memory can be replayed, and transferred to cortical regions for long-term memory storage [4][5][6]. During sleep spindles, which occur endogenously in temporal association with both the SO and

Cortical intrinsic currents
The pyramidal and interneuron cells of the cortex are modeled as two separate compartments (dendritic and axosomatic compartment) as initially proposed by Mainen and Sejnowski [25] based on Hodgkin-Huxley kinetics [24]. C m dV D dt ¼ À g L V D À E L ð Þ À g SD V D À V S ð Þ À I int D À I syn ; In Eq 1, C m and g L are the membrane capacitance and leakage conductance of the dendritic compartment. E L is the reversal potential, V D the dendric and V S is the axosomatic compartment membrane potential. g SD and g DS are the conductances between the axosomatic and dendritic compartments, respectively and g_SD = 1/(R � S soma � 165) and g_DS = 1/(R � S soma ) where R = 10 MO and S soma = 1.0 � 10 −6 cm 2 . I int D is the sum of active dendritic, I int S the sum of active axosomatic currents and I syn the sum of synaptic currents. I int D and I int S are the sum of the following intrinsic currents: I int S ¼ I Na þ I Na p ð Þ þ I K : Thereby, I Na represents the fast sodium current, I Na(p) the persistent sodium current, I LK the potassium leak current, I HAV the high-threshold calcium current, I Kca the slow calcium-   PLOS ONE dependent potassium current, I KM the slow voltage-dependent non-inactivating potassium current and I K represents the delayed rectifier potassium current. The IN cell compartments have the same intrinsic currents except for I Na(p) that is only included in PY cells, The ratio of dendritic area to somatic area was set to ρ = 165 in PY cells, and to ρ = 50 in IN cells. All voltage-dependent currents I c were simulated in the same fashion: where g c is the maximum conductance, m is an activation gating variable, M is the number of activation gates, h is an inactivation gating variable, and N is the number of inactivation gates. V is the corresponding compartment voltage and E c is the reversal potential. The dynamics of all gating variables were solved with the same equations: where x is a gating variable, x = m or h, the temperature-related term, Q T = Q ((T-32)/10) = 2.9529 where Q = 2.3, T = 36˚C. α x and β x are voltage-dependent transition rates. All individual intrinsic currents are described in Table 1 and their units and parametric values are described in Table 2.

Thalamic intrinsic currents
Two separate sub-thalamic networks for fast and slow spindles, respectively, were developed, each with a thalamocortical/relay (TC) and reticular (RE) cell layer. Cells of each layer were modeled based on a single compartment (somatic compartment) using the same voltagedependent and calcium-dependent currents dynamics as expressed by Hodgkin-Huxley kinetics schemes [24]: where C m is the membrane capacitance, g L the leakage conductance, E L the reversal potential, and V the voltage of the compartment. I syn denotes the sum of the synaptic currents and similarly I int the sum of the active intrinsic currents. The sum of these active currents for TC, I int TC and RE, I int RE are described as here I Na represents the fast sodium current, I K the fast potassium current [26], I KL the potassium leak current, I h the hyperpolarization-activated cation current [27], I T the low-threshold calcium current in TC [28] and I T in RE neuron [29]. The potassium leak current is I KL = g KL (V-E KL ) in both TC and RE cells where g KL is potassium leak conductance and E KL is the potassium reversal potential (E KL = -95 mV). Calcium dynamics for thalamic cells is described by; where [Ca] 1 = 2.4 � 10 −4 mM, A = 5.1819 � 10 −5 mM cm 2 /(ms μA) and τ = 5ms. All individual voltage-dependent currents were simulated in the same fashion as the cortical intrinsic currents given in Table 1. Table 2 gives the parametric values developed in our model.

Synaptic currents
For synaptic signaling four types of synaptic currents, I syn were used, three (AMPARs, GABA A Rs, and NMDARs) were modeled by the first-ordered activation scheme [30,31]. Accordingly, these synaptic currents are given by where g syn is the maximal synaptic conductance, where t 0 is the time for receptor activation and θ(x) is the Heaviside function [32]. The duration and amplitude parameters for the neurotransmitter pulse are t max = 0.03 ms and A = 0.5. The synaptic current rate constants for AMPARs were α = 1.1 ms and β = 0.19 ms, for GABA A Rs α = 10.5 ms and β = 0.166 ms, and for NMDARs α = 1 ms and β = 0.0067 ms. Intracortical currents were modified by multiplying the short-term depression term "D" [33,34] with the maximal synaptic conductance in Eq 3 for AMPARs and GABA A Rs receptors: where D is the amount of available synaptic resources, calculated by the following scheme: where the synaptic resources time recovery is τ = 700 ms, the interval between nth and (n+1) Δt, and the fraction of resources used for each action potential is U (for AMPARs U = .07, for GABA A Rs U = .073).
The fourth synaptic current, GABA B Rs is computed by a higher-ordered activation scheme that involves potassium channel activation by a G-protein, [30,35]: where [G] reflects the G-protein concentration, [R] the fraction of activated receptors, and the potassium reversal potential E K = -95 mV. K1 = 0.052 m M -1 ms -1 , K2 = 0.0013 ms -1 , K3 = 0.098 ms -1 , and k4 = 100μ M 4 were the rate constants. The maximal synaptic conductance used here for each synapse is described in Table 3.

Network geometry
The network model is comprised of six one-dimensional layers of neurons (Fig 1). Each layer of cells has N neurons, (N = 40) except the PY neurons layer, which has 5N neurons (200 neurons) [36]. The first and second cortical layer of PY and IN neurons initiate SOs. The third and fourth are thalamic layers for fast spindle initiation, and similarly the fifth and sixth are also thalamic layers that initiate slow spindles. The radii of synaptic connections between different layers are described in Table 3. For each SO cycle initiation, EPSPs and IPSPs miniature currents were implemented to PY-PY, PY-IN and IN-PY cells via AMPARs and GABA A Rs receptors [37]. These mini currents emerge~100 ms after the start of the SO downstate. For Poisson input implementation, NetStim.noise was set to 1 in the NEURON simulator.

Computational environment
The model was simulated in the NEURON 7.6 simulation environment [38] and it was run on a MacBook Pro 2015. For data analysis, MATLAB (R2020a) and eeglab tool were used.

Results
In this study, we developed a conductance based thalamocortical model for NREM sleep. As far as we know this is the first model for normal NREM sleep including fast and slow spindles written in the NEURON simulation environment. Our model exhibits both fast and slow spindles along with SOs. Both fast and slow spindles were initiated in two separate thalamic subnetworks. The main network comprises six layers of cells (Fig 1), the top two of which are cortical layers of PY and IN cells for SOs, while the middle two are thalamic layers of TC and RE cells that initiate the fast spindles ('fast thalamic subnetwork') and the lower two layers are for slow spindles ('slow thalamic subnetwork'). The thalamic subnetwork for slow spindles has the same type of intrinsic and synaptic current dynamics and network architecture as the fast thalamic subnetwork, however, it has a larger hyperpolarization than the fast thalamic subnetwork. Hyperpolarization of the slow thalamic subnetwork was increased by setting the reversal potential of GABA A Rs to E GABAA = -88 mV in TC cells. The reversal potential of passive currents was set to E L = -77 mV in TC cells and E L = -82 mV in RE cells.

Model initiation of SO and spindle generation
The main thalamocortical network is initiated by mini synaptic current to both cortical layers during the hyperpolarized Down state. Once the SO cycle is initiated, the mini synaptic current is terminated. During the Down state, the mini synaptic current activates the persistent sodium current of PY neurons and consequently these PY neurons depolarize and reach firing threshold. Initially only one or few PY neurons generate an action potential, yet they target their neighboring PYs by strong PY-PY excitatory connections. Due to the strong PY-PY excitatory connections and persistent sodium current they sustain this depolarized Up state for 500-1000 ms. The calcium dependent potassium current and progressive synaptic depression terminate the depolarized Up state and bring cortical network back to the Down state. After 100 ms of terminating the SO cycle, the mini synaptic current is again activated for the next SO cycle initiation and similarly this process is repeated for each SO cycle. Here we discuss the sequential flow of our network model after the activation of the cortical network. As the cortical network is activated, both layers of the fast thalamic subnetwork (TC and RE cells layers) also receive cortical inputs and become active. In the fast thalamic network, the interaction of TC and RE cells produces fast spindle oscillations with a major contribution of the TC hyperpolarization current (I h ) and transient calcium current (I T ) [39]. The fast thalamic output is sent back to the cortical network by TC cells. The cortical network receives this fast thalamic feedback~200 ms after the initiation of the SO cycle (Fig 2). The thalamic network for slow spindles receives cortical inputs with a 600 ms delay. We specifically chose this delay of 600 ms after the initiation of SO to achieve the generation of slow spindles during the second SO-half. The netcon class of the NEURON simulation environment was used to model the delay. Hyperpolarization of the slow thalamic network was increased by setting the reversal potential of GABA A Rs to E GABAA = -88 mV in TC cells. The reversal potential of passive currents was set to E L = -77 mV in TC cells and E L = -82 mV in RE cells. This more hyperpolarized thalamic network produced slow spindle oscillations that send excitation back to the cortical network via TC cells. The cortical network receives this slow thalamic input~800 ms after the initiation of the SO cycle (Fig 2 bottom two raster plots of entire slow thalamic network activity, left, and unit activity, right). Cortical LFP was calculated as the sum of the presynaptic currents

PLOS ONE
(AMPARs, NMDARs, GABA A Rs) of PY cells. The resultant cortical local field potential (LFP) is depicted in Fig 3A. The LFP was smoothed using the function; filter() (MATLAB) with a window size of 200 for better visualization of spindles, nested with the cortical LFP (Fig 3B). The fast thalamic input (Fig 3C and 3D, red) projected to the cortical network~200 ms after the initiation of the SO cycle, whereas the slow thalamic input (Fig 3C and 3D, green) projected to the cortical net-work~800 ms after the initiation of the SO cycle, nearly a few hundred ms before the termination of the SO cycle. Similarly, the fast thalamic inputs were found in every SO cycle but the weaker inputs did not have an impact on the cortical LFP (Fig 3D 2 nd to 4 th SO cycle).
Moreover, in the temporal window the thalamic inputs with a lifetime of less than 0.5 seconds were not considered proper spindles. In normal simulations, the number of fast spindles was high compared to that of slow spindles (see Fig 3C above). In our slow spindles results, generally the waning phase of the slow spindle was completed before the completion of the SO cycle, although sometimes a few spikes of the waning phase were also observed after the completion of the SO cycle, i.e. in the down state of the SO. Sometimes slow spindles initiated the second SO cycle before the completion of the preceding SO cycle (Fig 4B 6 th and 8 th SO cycle). Slow spindles initiated approximately 10% of SO cycles during normal or control simulation (See Fig 5I below).
For further investigation, mini synaptic currents that initiate SOs, were reduced to observe the role of the slow spindle in SOs initiation. The time period of SO down states was shorter in the presence of slow spindles (see Fig 5G-5I). The finding that slow spindles may contribute to maintenance of SO activity is the major finding of our model.
Independence of fast and slow spindle networks was investigated in another set of simulations in which network activity was produced after blocking the fast thalamic network (blocking layer 3 and 4), slow thalamic network (blocking layer 5 and 6) or by blocking both thalamic networks (blocking all four thalamic layers). Results indeed show the temporal properties of fast and slow spindles can be retained independently of one another (Fig 4). The thalamic layers were blocked by setting the synaptic conductance to zero between the cortical PY cells and thalamic cells (PY to TC and PY to RE cells).

Properties of SO-Spindle interactions
In-vivo fast spindles occur normally during the SO down-to-up transition and SO up state. We refer to this interval as the SO first-half in our model. Correspondingly, the second SOhalf (~800 ms after the initiation of the SO) characterizes the up-to-down state transition during which slow spindles occur. The fast thalamic subnetwork input projected to the cortical network increases the firing rate of cortical cells (Fig 5A and 5B). Their spiking occurred earlier as compared to spikes from non-thalamically innervated cortical cells (Fig 5D). Furthermore, the first SO-half associated with fast spindles obtained a larger amplitude (Fig 5E) than without input from the fast spindle thalamic network. The average amplitude of SO cycles coupled with fast spindles was~8.1 mV, whereas the average amplitude without spindles was~7.6 mV. Thus, in our model, fast spindles increase both the firing rate of cortical cells and SO amplitude. Slow spindles emerge later, during the second SO-half around~700 ms after the initiation of the SO cycle, when the SO has already reached its peak amplitude and starts to decline. Omission of the slow spindle thalamic subnetwork (Fig 5C vs. 5D) had two effects: a slowing of up state cortical firing rate, and a reduced SO duration. On the other hand, the duration of the second SO-half was more frequently longer in the presence of slow spindles than when only fast spindles were present (Fig 5C vs. 5A). Thus, in our model slow spindles appear to essentially assist in initiating the SO. To underpin this finding, we observed the

PLOS ONE
model output across a time span of 120 second. The total time spent in SO down states was only~16 seconds when slow spindles coupled to the SO up-to-down transition in the cortical LFP whereas without slow spindles a longer total time was spent in the down state, i.e.~32.6 seconds. In this control simulation approximately 10% of SO cycles were initiated by slow spindles (e.g., Fig 5I red circled SOs): In cases of SO initiation by slow spindles miniature synaptic currents were not required for the initiation of SO.
In the second-SO-half, the cortical network unit frequency is comprehensively reduced and cortical cells cease firing. To better distinguish between the interaction of slow and fast spindle subnetworks with SO properties in this critical second SO-half, we replaced slow spindles with the fast spindles (see Fig 6), i.e., fast spindles were generated both during the first and second SO-halves (Fig 6B). Fast spindles in the second SO-half eliminated the down/hyperpolarized state of SO, the cortical network remained in its state of increased firing (similar to the first SO-half), albeit exhibited activity was irregular compared to controlled simulation results ( Fig  6A). In another control simulation, the model only allowed for generation of fast spindles during the second SO-half. The down states of SOs again disappeared and inconsistent cortical LFPs were observed (not shown in figure). From these results, we conclude that slow spindles maintain a slow oscillation rhythm by preserving the SO down state. . PY firing rate is relatively increased by fast spindle input (red); C, with only slow spindles (fast spindles were blocked). PY firing rate is relatively increase by slow spindle input (green). D, without spindle input (both fast and slow spindles were blocked). E, Average cortical LFP amplitude in the presence of both fast and slow spindles (blue), fast spindles only (red), slow spindles only (green), and without spindles (gray). Average SO amplitude is lowest when spindles are absent. Asterisks show significant difference between the amplitudes (two sample t-test with � P < .01). F, Average SO cycle duration dependent upon presence of spindles. Duration was longest for the presence of slow spindles. Error bars describe standard deviation between SO cycles. Asterisks show significant difference between the SO cycle duration ( � P < .01). G, Resultant cortical LFP when the input of miniature synaptic current was reduced, but slow spindles were still generated by the model. The horizontal red bars indicate periods of increased down state duration due to weak mini synaptic current. H, Same as for G, however slow spindles are blocked. The down state durations are on average longer than when slow spindles are present (G). I, The cortical LFP with the input of miniature synaptic current (red) and slow spindles (green). The red circles indicate SO initiated by slow spindles without mini current. https://doi.org/10.1371/journal.pone.0277772.g005

PLOS ONE
The thalamic membrane potential level (in both TC and RE cell layers) differs, physiologically, between depths of NREM sleep and thus plays a role in the likelihood of spindle generation [40,41]. To investigate the response of the model, we increased the hyperpolarization level of fast thalamic subnetwork. Hyperpolarization was increased by altering the conductance of potassium leaked current (g KL ): g KL = .033 mS/cm 2 (from g KL = .03 mS/cm 2 ). The hyperpolarization response is depicted by raster plots of the all TC cells in the fast (40 TC cells) and slow (40 TC cells) thalamic subnetworks and by the power spectra of the corresponding cortical LFPs (Fig 7): The membrane potential of TC cells of the fast thalamic subnetwork decreased with increased thalamic hyperpolarization compared to control fast thalamic subnetwork activity (Fig 7A and 7B). In each fast thalamic subnetwork with hyperpolarization a decreased number of spikes per event (~6 spikes) were observed relative to control simulations (~9 spikes/event; Fig 7D). Fast spindle LFP power also decreased with hyperpolarization (Fig 7B, left). Interestingly, this increase in hyperpolarization of the fast thalamic subnetwork resulted in a stronger slow thalamic subnetwork input to the cortical network and subsequently, the number of slow thalamic subnetwork events was increased (~25 events per minutes) as compared to the control simulations (~19 events per minute: Fig 7C). The resultant cortical LFP power in the slow frequency range remained the same after hyperpolarization (Fig 7ab right). In summary, the hyperpolarization of the fast thalamic subnetwork reduces fast spindle activity and increases activity within the slow thalamic subnetwork.
Finally, we disintegrated all three subnetworks; cortical, fast thalamic, and slow thalamic subnetworks to validate the main rhythms; SO, fast spindles and slow spindles in each individual network. The cortical network initiation procedure was just like in the integrated model. Both fast and slow thalamic subnetworks were initiated by injecting step current to TC and RE cells. In both thalamic subnetworks, RE cells got .09 nA current whereas TC cells got .065 nA current for 600 ms after every 3 seconds (see Fig 8). All three networks successfully generated expected frequency range. Cortical network generated LFP with little high frequency~1.3 Hz (Fig 8A). Thalamic results were observed in single cell activity. TC cells of the fast thalamic subnetwork fires with frequency of~16 Hz. Moreover, the fast thalamic subnetwork shows  weaker response with high potassium leak current (~14 Hz) (Fig 8B and 8C). Similarly, the slow thalamic subnetwork cells fire with frequency~10 Hz (Fig 8D).
In summary, in our model slow spindle properties can be modeled by a thalamic subnetwork that differs from the fast spindle thalamic subnetwork; Our model revealed the following properties of SO and spindle interactions: 1. fast spindles increased the SO amplitude; 2. slow spindles assisted in initiating the SO; 3. slow spindles regularized and stabilized the SO rhythm by preserving the SO down state; and 4. hyperpolarization of the fast thalamic subnetwork reduced fast spindle density and increased the number of slow thalamic network events (slow spindles).

Discussion
Sleep spindles are the other prominent EEG rhythm of NREM sleep. These are initiated by the interaction of thalamocortical (TC) and reticular (RE) cells of thalamic nuclei [39]. These

PLOS ONE
spindles are also sent to the cortical network by collaterals of TC cells. Generally, these are nested with cortical SOs at the initial phase of the SO cycle, the first-half phase in primates and rodents. In human sleep studies, spindles are classified into two categories, namely fast (12-16 Hz) and slow (8)(9)(10)(11)(12) spindles. The fast spindles are mostly found in central and parietal cortical regions, and these are also nested with SOs during the first-half of up state. Slow spindles are observed in the frontal cortical area, and these are nested with SO during the secondhalf of up state, [18,42]. Fast spindles are extensively studied in the experimental and computational arena regarding their role in long-term potentiation and synaptic plasticity [8]. Conversely, slow spindles need lot of work to find clear answers about their origin, cellular mechanism, and most importantly their role in such a framework (SO, fast spindles, hippocampal ripples) that is dedicated for memory consolidation. As a step ahead, it would be useful to explore the functional relation between SOs and slow spindles.
In this study, we developed a thalamocortical model for SOs, fast and slow sleep spindles. The novel aspect is the inclusion of thalamic-based slow spindle generation. Two independent sub-thalamic networks were developed for each fast and slow spindles. Physiologically fast and slow spindles occur at the SO down-to-up transition / SO up-phase and at the SO up-to-down transition, respectively. In the model fast spindles are nested with SOs during the first-half of the SOs cycle (SO down-to-up transition / SO up-phase) and the slow spindles are nested with SOs during the second SO-half (end of the SO up state/up-to-down transition). Slow spindles are produced with the same intrinsic and synaptic current dynamics as fast spindles but the thalamic subnetwork is in a more hyperpolarized state. According to our model results, slow spindles can initiate SOs and may thereby facilitate the maintenance of ongoing SOs. A contribution of slow spindles to the maintenance of activity in the SO state could also be deduced from the decrease in SO duration on omission of slow spindles. Faster electrophysiological activity ensues as NREM sleep lightens.
Fast spindles could not replicate these actions of slow spindles on the SO. When simulations were run for fast spindles nesting within the second SO-half (end of the up-state), the cortical network exhibited an irregular LFP pattern and the down states of SOs practically disappeared. Moreover, the fast thalamic subnetwork response was reduced when the hyperpolarization level in this network was increased, in contrast to the facilitatory response of the slow thalamic subnetwork to hyperpolarization. Thus, our model shows that slow spindles can facilitate cortical network activity while maintaining the natural rhythm of SOs.
To which neurophysiological processes are the model properties consistent? Prominent initial studies revealed that SO arise from layer 5 pyramidal cells and that SOs continue despite thalamic deafferentiation [43] e.g., through intrinsic activity emerging in layer 5 pyramidal cells [11,43,44]. It was recognized that thalamic input can contribute to cortical Up-states [45], however, the argument that SOs do not require thalamic input for initiation contrasts more recent studies that disclosed thalamic activity preceding the onset of cortical Up states, and also that severing thalamocortical connections reduced the incidence of spontaneous cortical Up states [46], and also modified ongoing SO frequency [47]. Experiments on anaesthetized cats demonstrated that thalamic oscillations contribute importantly to the cortical network in generating SOs, and can results explain previous contradictory findings [48], reviewed in [49]. Thalamic output at the time of slow spindle activity may thus deliver the input to the cortex required for triggering cortical Up states. Some thalamocortical cells, including the ventral lateral posterior (VLp) nucleus as a possible thalamic source of slow frequency spindles [22], possess both core-like and matrix neurons, and project thus to both superficial and deeper cortical layers, including axonal arborizations to layer V [50][51][52]. The intrinsic initiation of SOs in cortical layer 5 is discussed in the context of one the of two possible mechanisms, firstly, by persistently active pacemaker-like cortical cells, and secondly by temporal summation of spontaneous synaptic activity [44,49]. Any one of both mechanisms cannot initiate the SOs unless it counters the activity-dependent K + conductances that are activated during active states. Consistent with our modeling results slow spindle could assist in SO initiation by providing some depolarization to counter activity-dependent hyperpolarization conductances. Independent on whether such bilayer cortical input could be beneficial for SO initiation in layer 5 pyramidal cells, the above cytoarchitectonics and neurophysiological processes may allow to explain model results on slow spindle-SO initiation.
Processes underlying the maintenance of the SO rhythm are undoubtedly even more diverse and complex than can be reflected in our model. Human intracranial recordings and rodent local field potentials describe mechanisms on how thalamic spindles may drive cortical spindles emerging during the SO down-to-up transition [53,54]. Slow and fast spindles are found however to occur during the same SO phase [53], thus findings on the temporal relationship between thalamic and cortical down states are not comparable with our modeling results. These inconsistencies in slow spindle timing present a major focus of ongoing research [55][56][57].
Modeling results presented two properties regarding the impact of fast spindle on the slow spindle network. Firstly, that fast spindles increased the firing rate of cortical cells and SO amplitude, and secondly that the hyperpolarization of the fast thalamic network reduced fast spindle, yet increased slow spindle activity. Both properties are associated with membrane potential level of thalamic cells, and coincide with results linking increased thalamic hyperpolarization to the emergence of delta oscillations and increased NREM sleep depth [58]. In deeper NREM sleep occurrence and frequency of spindles are also reduced in humans [21].
Taken together, our modeling results present above most a further platform to test the potential role to thalamic spindles in SO initiation. In future work, we will extend this study by making structural and physiological modifications for more detailed experiments. For instance, by including neuronal plasticity in the current model we could analzye the potential contribution of slow spindles to this function.